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ABSTRACT 

This  paper  summarizes  a  new  capability  in  the  topology  optimization,  which  employs  an 
advanced  substructuring  technique  in  the  multi-domain  topology  optimization  for  the  purpose  of 
designing  a  complex  and  practical  mechanical  system.  In  the  substructuring  method,  a  structural 
system  is  divided  into  a  series  of  substructures,  which  can  be  designed  separately  with  the 
coupling  of  the  subsystems.  The  condensation  technique  is  applied  to  the  substructures  where 
the  design  changes  are  not  considered  at  the  current  step.  For  the  substructures  where  the 
optimal  design  is  sought,  the  multi-domain  topology  optimization  is  applied  in  order  to  achieve 
the  design  objectives  prescribed.  The  embedded  advanced  substructuring  technique  allows  users 
to  select  nodal  displacements  and  modal  coordinates  as  generalized  coordinates  in  the  reduced- 
order  models,  and  to  use  advanced  Quasi-Static  modes  for  high  efficient  modal  analyses. 
Examples  are  illustrated  in  order  to  demonstrate  the  feasibility  of  the  proposed  technique  to  the 
topology  optimization  problem. 
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1.  INTRODUCTION 

Topology  optimization  (1-4)  has  been  widely  adopted  in  the  design  process  for  developing  light¬ 
weight  and  high-performance  structures,  especially,  in  automotive  vehicle  structural  systems.  A 
bottleneck,  however,  is  how  to  apply  the  topology  optimization  technique  to  deal  with  a  large 
and  complex  engineering  structure,  such  as  a  major  subsystem  in  a  full  vehicle  system.  In  the 
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practical  design,  a  structural  system  may  contain  a  large  number  of  components  and 
substructures  and  each  substructure  may  have  a  large  number  (hundreds  of  thousands)  degrees  of 
freedom  in  its  discretized  finite  element  model.  The  direct  adoption  of  the  whole  system  in  a 
topology  optimization  process  is  inefficient  and  even  formidable  from  the  aspect  of 
computational  resource  and  the  aspect  to  find  a  solution  for  the  optimization  problem. 

This  research  aims  at  the  implementation  of  advanced  substructuring  techniques  (11,16-19)  in 
the  multi-domain  topology  optimization  for  designing  complex  mechanical  systems.  In  the 
substructuring  method  (17-19),  a  mechanical  system  is  divided  into  a  series  of,  or  multileveled, 
substructures,  which  can  be  designed  separately  with  the  coupling  effects  of  the  subsystems.  The 
advanced  substructuring  technique  allows  users  to  select  nodal  displacements  at  the  interfaces  of 
the  substructures  and  modal  coordinate  as  generalized  coordinates  in  the  reduced-order  model, 
and  to  use  advanced  Quasi-Static  modes  (11,17)  for  high  efficient  modal  analyses.  The 
condensation  technique  can  be  applied  to  the  substructures  where  the  design  change  is  not 
considered  at  the  current  step.  For  the  substructures  where  the  optimal  design  is  sought,  the 
multi-domain  topology  optimization  (20)  can  be  applied  for  those  subdomains  in  order  to 
achieve  the  design  objectives  of  the  whole  structural  system.  Another  usage  of  this  approach  is 
that,  the  topology  optimization  can  be  zoomed  into  any  subdomain  within  a  design  domain  at 
anytime  in  order  to  lay  out  the  details  of  the  design  at  local  area,  or  to  reduce  the  stress 
concentration  at  a  localized  location.  Using  this  divide-and-conquer  design  methodology,  an 
optimal  structural  design  can  be  obtained  for  the  coupled  substructures  with  reduced 
computational  cost  and  better  solvability  of  the  design  problem  of  large  size,  which  are  crucial  in 
the  application  of  the  topology  optimization  technique  for  real  engineering  design  problems. 

2.  MULTI-DOMAIN  TOPOLOGY  OPTIMIZATION 

A  multi-domain  topology  optimization  (MDTO)  approach  is  proposed  in  (20),  which  can  be 
considered  as  a  generalization  of  the  standard,  single-domain  topology  optimization  developed  in 
the  earlier  stage  (5-9).  In  contrast  to  single-domain  topology  optimization,  in  which  a  given 
amount  of  the  material  is  assigned  to  the  entire  design  domain,  MDTO  allows  the  designer  to 
assign  different  amounts  of  the  material,  or  even  different  materials,  to  the  various  subdomains 
of  the  structure.  Figure  1  depicts  a  structural  domain  that  is  divided  into  three  subdomains,  where 
a  certain  amount  of  material  A  is  assigned  to  Subdomain  1 ,  a  different  amount  of  material  B  is 
distributed  into  Subdomain  2,  and  Subdomain  3  is  considered  as  a  non-design  domain,  for  which 
the  material  distribution  is  not  allowed  to  change  at  the  current  design  stage. 


Figure  1.  A  multi-domain  topology  optimization  problem 
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In  the  general  case,  the  MDTO  problem  can  be  formulated  as 
Minimize  / (X) 

Subject  to:  h  (X)  < 0  (J  =  l,2,...,m) 

(1) 

Xj  <xt<Xi  (?  =  1,2, ...,«) 

(and  State  Equations) 

where  f=J[X)  denotes  the  objective  function,  A.  =  A;(X)  denotes  the  y-  th  constraint  function  for 
the  volume  (or  weight)  of  the y'-th  substructure  in  the y'-th  subdomain  (where  y  =  1,  2,  •  •  •,  m  );  m  is 
the  total  number  of  the  subdomains;  X  =  {x1?  x2,  ...,  xn)T  denotes  the  vector  of  the  design 

variables,  where  n  is  the  total  number  of  the  design  variables;  and  xt  and  x,  are  the  lower  and 
upper  bounds  of  design  variable  xn  respectively.  Note  that  f[X)  in  Eq.  (1)  also  needs  to  satisfy 
the  state  equations  for  the  structural  problem  at  hand.  These  state  equations  may  include,  for 
example,  the  static  equilibrium  equation,  the  equation  that  defines  the  free  vibration  eigenvalue 
problem,  or  the  equation  for  the  dynamic  forced  response. 

The  basic  idea  in  the  topology  optimization  technique  proposed  by  Bendsoe  and  Kikuchi  in  1988 
(4)  is  to  transform  the  optimal  topology  design  problem  into  an  equivalent  Optimal  Material 
Distribution  (OMD)  problem,  as  illustrated  in  Fig.  2.  Here  the  structural  domain  is  assumed  to  be 
filled  with  a  non-homogeneous  composite  material  characterized  by  a  variable  microstructure.  A 
typical  microstructure  is  formed  inside  an  empty  rectangle  in  a  unit  cell  with  three  design 
variables  a,  b  and  0,  which  are  dimensions  and  orientation  of  the  microstructure,  respectively. 
Using  homogenization,  the  effective  material  constants  can  be  obtained  as  (4): 


D"  =D H  (a,  b,  0 ) ,  pH  =  pH  (a,  b )  (2) 

where  D"  denotes  the  effective  elastic  coefficient  matrix  at  the  material  point,  and  p"  is  the 
associated  effective  mass  density.  Both  D"  and  pH  are  functions  of  the  design  variables  a  and  b, 
and  D"  is  also  a  function  of  the  orientation  variable  Oat  the  material  point. 

In  the  optimization  process,  the  microstructure  can  vary  anywhere  between  "empty"  (a=b= 0)  and 
"solid"  (a=b=\ )  using  the  design  variables  a  and  b,  and  it  can  be  rotated  using  the  orientation 
variable  0.  Therefore,  if  one  assumes  that  the  total  amount  of  the  material,  which  is  prescribed 
for  the  design  problem  at  hand,  remains  constant  in  the  optimization  process,  then  the  material 
will  be  moved  from  a  region  of  the  structural  domain  into  another  region  to  produce  a  new 
distribution  of  the  material.  Structural  optimization  technique  (10,12-15)  has  been  used  to  search 
for  the  optimal  material  distribution.  By  moving  and  orienting  the  material  so  as  to  improve  the 
objective  function  of  the  optimization  problem,  one  can  finally  obtain  an  OMD  that  corresponds 
to  the  optimal  structure.  Various  microstructure  models  have  been  developed  based  on  this 
fundamental  concept,  including  artificial  material  models,  which  represent  a  simplification  of 
this  concept  (1-4). 
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Unit  cell 


Figure  2.  Basic  concept  of  the  homogenization-based  topology  optimization  method 

3.  ADVANCED  SUBSTRUCTURING  METHODS 

Using  a  standard  finite  element  method,  the  state  equation  governing  the  dynamic  response  of 
the  z-th  structure  shown  in  Fig.  1  can  be  obtained  as: 

M(i)u(0  +  C(,)u(i)  +  K(i)u(0  =  f (0  (3) 

where  M(<),  C(l\  and  K(,)  are  the  mass,  damping,  and  stiffness  matrices  of  the  z'-th  substructure  in 
the  system,  respectively.  Mw  and  Kw  can  be  obtained  by  assembling  the  elementary  mass  and 
stiffness  matrices  as 


nel  nd 

M(0  =  me  and  K<0  =  ke  (4) 

e=l  e=l 

where  m  =m  (a.  b)  and  k  =k  (a,  b  ,  6)  are  the  mass  and  stiffness  matrices  of  the  finite 
element  e,  which  are  functions  of  ae,  be,  and  de ,  where  ae,  be,  and  0e  are  the  discretized  design 
variables  of  a,  b,  and  6  ,  respectively,  at  the  elementary  level  (e  =  1,  2,  ■■■,  nel).  In  the  special 
case  of  proportional  damping,  the  viscous  damping  matrix  C(/>  in  Eq.  (3)  can  be  obtained  as 
C(,)  =  /?,K(,) ,  where  /?  is  the  damping  coefficient  of  the  z'-th  substructure. 

State  equations  for  the  topology  optimization  problem,  Eq.  (3),  usually  feature  a  very  large 
number  of  variables  (typically  the  nodal  displacements).  Even  for  a  simple  structure,  a  finely 
discretized  finite  element  model  can  involve  from  thousands  to  millions  of  nodal  variables.  In 
topology  optimization,  such  fine  meshes  are  usually  required  for  obtaining  smooth  boundaries 
and  interfaces.  Substructuring  methods  provide  an  effective  tool  for  condensing  the  analysis 
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variables  into  a  much  smaller  set,  thus  greatly  reducing  computer  memory  requirements  and 
increase  computational  efficiency.  This  approach  is  particularly  natural  and  useful  for  multi- 
domain  topology  optimization  problems.  While  substructuring  can  be  used  in  several  ways  to 
improve  computational  efficiency,  in  this  paper  the  condensed  analysis  variables  are  those 
associated  with  the  non-design  subdomains  of  the  structure  (see  Fig.  1).  The  nodal  displacement 
variables  can  be  condensed  into  those  at  the  boundaries  of  the  subdomain  when  a  static  analysis 
is  considered,  and  with  an  additional  small  number  of  generalized  modal  coordinates  when  a 
dynamic  problem  is  considered.  Note  that  the  design  and  non-design  subdomains  can  be 
switched  if  necessary  during  the  different  design  stages. 

In  the  case  of  a  static  condensation,  assuming  that  ua  denotes  the  vector  of  nodal  displacements 
associated  with  the  active  nodes  (usually  the  boundary  nodes)  of  a  non-design  subdomain,  and 
uo  denotes  the  vector  of  the  nodal  displacements  associated  with  the  other  nodes  (usually  the 
internal  nodes)  of  the  same  subdomain,  the  state  equation  for  the  subdomain  can  be  written  as 


where  k..(7,  j  =  o,  a)  are  the  blocks  of  the  stiffness  matrix  of  the  subdomain  associated  with  u0 
and  ua ,  respectively,  and  f  and  f,  are  the  corresponding  nodal  force  vectors.  Then  Eq.  (5)  can 
be  condensed  into  an  equation  that  involves  only  ua ,  namely, 

Ka*a=Z  (6) 

where  k*  =  k  -k  k“'k  and  f*  =  f  -k  k“'f  .  For  the  sake  of  simplicity,  but  not  losing 

aa  aa  ao  oo  oa  a  a  ao  oo  o  r  J  ?  o 

the  generality,  in  the  following  discussions,  we  assume  fo  =  0 ,  and  therefore  we  have  fj  =  ia . 
Note  that  the  coordinate  transformation  for  the  condensation  is 


where  the  transformation  matrix  is 


(V) 


(8) 


Each  column  of  C  is  called  a  “ static  mode ”  or  a  “ constraint  mode ”. 

In  the  case  of  a  dynamic  condensation,  the  state  equation  can  be  written  as 
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where  m tJ(i,j  =  o,  a)  denote  the  blocks  of  the  mass  matrix.  We  now  divide  the  active  nodal  set 
(a-set)  into  two  subsets  p- set  and  5-set  (namely  a  =  p®s),  where  the  5-set  is  defined  by  those 
nodes  that  will  be  fixed  in  the  eigenmodes  extraction  process,  and  the  p- set  is  defined  by  the 
other  nodes  (those  will  be  free  in  the  eigenmodes  extraction  process).  Defining  an  eigenvalue 
problem  with  respect  to  the  nodal  set  n  =  o® p ,  namely 


rrk 

k 
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m 

M=o 

oo 

op 

-k 

oo 

op 
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K„J 

vL  p° 

pp  _ 

L  p° 

pp  J 

) 

(10) 


we  now  can  obtain  a  set  of  normal  modes  from  Eq.  (10),  namely  O  = 


O 

O 


on 

pn 


Using  the  normal  modes  calculated  from  Eq.  (10),  the  transformation  equation,  Eq.  (8)  is 
extended  as 


\Uo 

I  u„ 


D 


u„ 


(11) 


where  q„  is  a  vector  of  so-called  generalized  modal  coordinates.  The  transformation  matrix  D 
for  the  dynamic  problem  can  then  be  obtained  as 


D 


The  matrix  Y  = 


®  -Y  ®  Y 

on  op  pn  o 


(ll) 


Y„ 


contains  the  so-called  “ quasi-static  modes'”  (QSM)  (11,17), 


where  Y  =  [  Y0/,  Y0J  ] ,  which  can  be  calculated  by  solving  the  quasi-static  (frequency  response) 
problem  associated  with  Eq.  (2).  Theoretically,  Y  can  be  written  as 


Y  = 


-(k™  -  °Kmoo  )  1  (k0a  -  0Kmoa  ) 


(12) 


where  coc  is  the  “ central  frequency ,”  which  is  determined  by  the  frequency  range  of  interest 
(11,17).  Note  that  the  size  of  q„  in  Eq.  (10)  is  usually  far  smaller  than  u0 ,  therefore,  the 
coordinate  transformation  Eq.  (11)  (as  well  as  Eq.  (8))  significantly  reduces  the  size  of  the 
analysis  problem.  For  a  static  analysis,  Eq.  (8)  does  not  induce  any  error  into  the  original 
analysis  problem.  But  for  a  dynamic  analysis  problem,  Eq.  (11)  is  an  approximation,  and  the 
error  can  be  controlled  by  properly  choosing  the  central  frequency  and/or  considering  additional 
modal  coordinates  in  Eq.  (11).  Note  that  for  coc=  0,  the  QSM  defined  in  Eq.  (12)  reduce  to  the 
traditional  static  (constraint)  modes  proposed  by  Hurty  (18)  and  Craig  and  Bampton  (19),  which 
are  given  in  Eq.  (8).  Compared  with  the  static  modes,  the  QSM  can  significantly  reduce  both  the 
size  of  the  analysis  problem  and  the  error  induced  by  the  coordinate  reduction  process, 
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particularly  within  the  frequency  range  of  interest.  This  is  because  not  only  the  higher-frequency 
modes  are  truncated,  but  also  the  lower-frequency  modes  outside  the  frequency  range  of  interest. 
Therefore,  the  QSM  can  handle  dynamic  response  problems  within  higher  frequency  ranges,  and 
they  include  the  traditional  static  modes  as  a  special  case.  Note  that  the  central  frequency  is 
typically  selected  at  the  middle  of  the  frequency  domain  of  interest,  and  that  multiple  central 
frequencies  can  also  be  used  (17).  Finally,  note  that  if  p- set  defined  in  Eq.  (10)  is  null,  then  Eq. 
(11)  is  reduced  to  a  fixed-interface  CMS  method  and  the  Craig-Bampton  method  can  be 
considered  as  a  special  case  of  it  when  coc=0.  On  the  other  hand,  if  p-set  becomes  n-set,  the  Eq. 
(1 1)  is  reduced  to  a  so-called  free-interface  CMS  method. 

4.  EXAMPLES 

4.1  Refine  design  in  a  zoomed-in  area  A  cantilevered  sandwich  beam  design  problem  is  shown 
in  Fig.  3  to  illustrate  the  multi-domain  topology  optimization  along  with  the  substructuring 
technique.  The  external  force,  F,  represents  a  combination  load  of  pushing  in  axial  direction  and 
bending  in  lateral  direction.  Four  design  subdomains  are  specified  in  Fig.  3,  and  they  are  shown 
in  green  and  yellow.  The  blue  non-design  domains  are  filled  with  solid-cell  material  (a=b=  1  in 
Eq.  2).  The  design  objective  is  to  obtain  an  optimal  material  distribution  in  the  design  domain  so 
that  the  cantilevered  beam  has  the  maximum  rigidity  under  the  specified  external  load,  F. 


Figure  3.  A  cantilevered  sandwich  beam  design  problem. 

A  multi-domain  topology  optimization  problem  is  considered  with  a  finite  element  model  of 
3,906  degrees-of- freedom.  The  material  density  for  each  design  subdomain  is  given  as  20%  with 
an  even  distribution  among  the  subdomains. 

A  converged  optimal  material  distribution  for  the  first  design  problem  is  shown  in  Fig.  4-a.  The 
von  Mises  stress  resulted  from  this  design  is  shown  in  Fig.  4-a,  with  the  normalized  maximum 
stress  16.6,  occurring  at  the  elements  near  subdomain  1  and  2.  In  order  to  reduce  the  maximum 
stress  in  the  first  design,  additional  13%  material  is  added  to  subdomains  1  and  2  while  other 
design  domains  are  fixed  at  the  current  step.  The  substructuring  technique  is  utilized  to  condense 
the  DOFs  in  subdomains  3  and  4,  which  reduce  the  total  DOFs  from  3,906  to  960.  The  design 
domain  in  the  second  step  optimization  is  shown  in  the  grey  shadow  in  Fig.  4-b. 

After  the  second  step  topology  optimization,  the  maximum  von  Mises  stress  is  then  reduced  to 
13.7  as  shown  in  Fig.  4-c.  To  further  reduce  the  maximum  stress  in  the  current  design,  the 
topology  optimization  is  zoomed  into  two  new  small  design  domains  near  the  stress 
concentration  areas  as  shown  in  Fig.  4-d  in  grey  color.  All  the  remaining  subdomains,  the 
original  design  subdomains  1-4  and  the  blue  non-design  domain  in  Fig.  3,  are  condensed  through 
the  substructuring  technique.  The  total  DOFs  of  the  FE  model  for  the  final  stage  optimal  design 
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becomes  336,  which  is  significantly  smaller  than  the  original  3,906  in  the  original  optimization 
problem.  The  final  optimal  design  is  shown  in  Fig.  4-e.  The  maximum  von  Mises  stress  is  now 
reduced  to  10.4  from  the  maximum  stress  16.6  in  the  first  step  optimization  by  adding  a  small 
percentage  of  the  material. 


Maximum  stress  =  16.6 


(a)  von  Mises  stress  after  the  1st  topology  optimization 


(b)  design  problem  for  the  2nd  design  process 


(c)  von  Mises  stress  after  the  2nd  topology  optimization 


(d)  design  problem  for  the  3rd  design  process 


Maximum  stress  =  10.4 


(e)  von  Mises  stress  after  the  final  topology  optimization  (f)  the  final  optimum  design  with  smoothed  boundaries 

Figure  4.  The  optimization  design  procedure  and  the  results 

Through  the  three- stage  multi-domain  topology  optimization,  an  optimal  design  has  been 
obtained  with  the  maximum  rigidity  and  significantly  reduced  maximum  stress,  the  later 
improves  the  strength  of  the  structure.  The  maximum  stress  is  reduced  with  the  substructuring- 
technique-based  refine  design  procedure  developed  in  this  paper. 

4.2  Coupled  substructures  design  The  second  example  demonstrates  a  follow-up  optimal 
design  problem  for  the  cantilevered  beam  illustrated  in  section  4. 1 .  Given  the  optimal  design 
layout  shown  in  Fig.  4-f,  we  now  consider  an  optimal  fixture  design  problem  for  the  structure. 
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The  left-end  fixture  part  is  redesigned  using  the  topology  optimization  process  with  reduced 
amount  of  material  for  the  fixture.  The  initial  design  problem  is  shown  in  Fig.  5-a,  while  light 
blue  color  indicates  the  design  domain  for  the  fixture. 

In  Fig.  5-a,  the  material  constraint  for  fixture  design  domain  is  37.5%.  The  optimal  structure 
obtained  in  section  4.1  is  condensed  with  substructuring  technique.  The  final  DOFs  of  the  FE 
model  used  in  current  optimization  is  914,  while  the  DOFs  of  the  full  FE  model  is  3,906.  The 
design  result  for  the  fixture  is  shown  in  Fig.  5-b,c,d.  It  is  seen  that  materials  in  the  design 
domain  are  moved  to  the  upper  and  lower  edges  to  withstand  the  lateral  bending  load  and  a  hole 
is  opened  in  the  middle  in  order  to  reduce  the  weight  of  the  fixture.  The  final  optimal  design  is 
shown  in  Fig.  5-d  with  the  smoothed  boundaries. 


(a)  fixture  design  problem  with  the  fixed  structural  (b)  optimal  fixture  design  obtained  with  the  fixed 
design  structural  layout 


(c)  stress  distribution  and  displacement  of  the  (d)  final  shape  of  the  system  with  the  smoothed 

coupled  structure-fixture  system  boundaries 

Figure  5.  The  topology  optimization  in  the  fixture  region. 


5.  CONCLUSIONS 

An  optimal  design  problem  is  illustrated  to  demonstrate  the  benefits  from  embedding  the 
substructuring  technique  in  the  multi-domain  topology  optimization  process.  It  is  seen  that  the 
substructuring  technique  can  reduce  the  size  of  the  finite  element  model,  thus  improves  the 
solvability  of  the  design  problem.  Stress  concentration  can  be  overcome  efficiently  using  the 
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zoomed-in  topology  optimization  in  the  neighborhood  of  the  stress  concentration.  The  proposed 
substructuring  technique  can  be  also  applied  to  dynamic  problems  such  as  the  eigenmode 
optimization,  and  the  frequency  response  optimization,  and  be  readily  extended  to  deal  with 
three-dimension  design  problems  for  complex  mechanical  systems. 
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